#ifndef __tri_hh__
#define __tri_hh__

class Tri
{
public:
  Tri(const double t[3*3]) : vertices(t) { compute_volume() ; }
  virtual ~Tri() {}
  virtual double* evaluate_basis ( int n_points, const double points[] ) = 0 ;
  virtual double* evaluate_gradient ( int n_points, const double points[] ) = 0 ;
  virtual int test() = 0 ;
  
//  virtual double* reference_to_physical_coords ( int n_points, double ref[] ) ;
//  virtual double* physical_to_reference_coords ( int n_points, double phy[] ) ;
  
//  double* barycentric_coords ( double p[3] ) ;
  virtual double element_volume () ;
  
protected:
  const double* vertices ;
  double volume;

private:
  void compute_volume() ;

};

inline
void Tri::compute_volume()
{
  volume = vertices[0] * ( vertices[4] - vertices[7] ) 
       + vertices[3] * ( vertices[7] - vertices[1] ) 
       + vertices[6] * ( vertices[1] - vertices[4] );
  
  // check if vertices are degenerate
  if ( volume < 1e-8 )
  {
    cout << "\n";
    cout << "Triangle3 - Fatal error!\n";
    cout << "  Element has zero or negative area.\n";
    cout << "  Area = " << volume << "\n";
    cout << "\n";
    if (volume < 0)
    {
      cout << "  Reorientation of node(2) necessary to achieve positivity.\n";
      cout << "  Physical Nodes:\n";
      cout << "\n";
      for (int j = 0; j < 3; j++ )
      {
        cout << "  "
            << setw(10) << vertices[0+j*3] << "  "
            << setw(10) << vertices[1+j*3] << "  "
            << setw(10) << vertices[2+j*3] << "\n";
      }
      volume = -volume;
    }
    else
      exit ( 1 );
  }

}

inline
double Tri::element_volume()
{
  return this->volume ;
}

class Tri3 : public Tri
{
public:
  Tri3(const double t[3*3]) : Tri(t) {} 
  virtual double *evaluate_basis ( int n_points, const double points[] ) ;
  virtual double* evaluate_gradient ( int n_points, const double points[] ) ;
  virtual int test() ;
  
private:
};



class Tri4 : public Tri
{
public:
  Tri4(const double t[4*3]) : Tri(t) {} 
  virtual double *evaluate_basis ( int n_points, const double points[] ) ;
  virtual double* evaluate_gradient ( int n_points, const double points[] ) ;
  virtual int test() ;
  
private:
};


class Tri6 : public Tri
{
public:
  Tri6(const double t[4*3]) : Tri(t) {} 
  virtual double *evaluate_basis ( int n_points, const double points[] ) ;
  virtual double* evaluate_gradient ( int n_points, const double points[] ) ;
  virtual int test() ;
  
private:
};




#endif // #ifndef __tri_hh__
